rice (Rice / Rician distribution)#
The Rice (Rician) distribution is a continuous distribution on nonnegative real numbers. It arises as the distribution of the magnitude of a 2D Gaussian vector with a nonzero mean (equivalently, the amplitude of a complex Gaussian random variable with a deterministic “line-of-sight” component).
1) Title & Classification#
Item |
Value |
|---|---|
Name |
Rice ( |
Type |
Continuous |
Support |
\(r \in [0,\infty)\) |
Parameters |
noncentrality \(\nu\ge 0\), scale \(\sigma>0\) |
Useful reparam. |
\(K = \nu^2/(2\sigma^2)\) (Rician \(K\)-factor) |
What you’ll be able to do after this notebook#
connect Rice to a shifted 2D Gaussian and to Rayleigh / noncentral \(\chi\) distributions
write the PDF and CDF in standard special-function notation
compute mean/variance/skew/kurtosis from raw moments
derive the PDF (polar change of variables) and the log-likelihood
sample from Rice using NumPy only
use
scipy.stats.ricefor evaluation, simulation, andfit
import numpy as np
import plotly.graph_objects as go
import os
import plotly.io as pio
from scipy import stats, special
from scipy.optimize import minimize
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=6, suppress=True)
rng = np.random.default_rng(42)
2) Intuition & Motivation#
What this distribution models#
A canonical construction:
Let \((X,Y)\) be independent with $\(X \sim \mathcal{N}(\nu,\sigma^2),\qquad Y \sim \mathcal{N}(0,\sigma^2).\)$
Define the amplitude (radius) $\(R = \sqrt{X^2 + Y^2}.\)$
Then \(R\) has a Rice distribution with parameters \((\nu,\sigma)\).
Intuition:
When \(\nu=0\), the mean is at the origin, so the radius is the familiar Rayleigh distribution.
When \(\nu>0\), the Gaussian cloud is shifted away from the origin; radii concentrate around \(\approx \nu\).
The dimensionless ratio \(K = \nu^2/(2\sigma^2)\) (often called the Rician \(K\)-factor) measures how strong the deterministic component is relative to the random scatter.
Typical real-world use cases#
Wireless communications: Rician fading amplitudes (line-of-sight + multipath scatter).
Radar / sonar: return amplitude with a deterministic component plus noise.
MRI magnitude images: magnitude of complex Gaussian noise leads to Rician-like likelihoods.
Signal processing: envelope of a sinusoid in additive Gaussian noise.
Relations to other distributions#
Rayleigh: \(\nu=0\) gives Rayleigh\((\sigma)\).
Noncentral \(\chi\): Rice is the noncentral-\(\chi\) distribution with degrees of freedom \(k=2\).
Noncentral \(\chi^2\): \(R^2/\sigma^2 \sim \chi'^2(2,\lambda)\) with noncentrality \(\lambda=(\nu/\sigma)^2\).
Approx. Normal for large \(K\): when \(\nu\gg\sigma\), \(R\) is close to \(\mathcal{N}(\nu,\sigma^2)\) truncated at \(0\) (the truncation becomes negligible).
3) Formal Definition#
We use the parameterization \((\nu,\sigma)\) with \(\nu\ge 0\) and \(\sigma>0\).
PDF#
For \(r\ge 0\):
and \(f(r)=0\) for \(r<0\). Here \(I_0(\cdot)\) is the modified Bessel function of the first kind (order 0).
A useful integral identity (often used in derivations) is:
CDF#
A standard special-function expression uses the Marcum \(Q\)-function of order 1:
One definition of \(Q_1(a,b)\) is
SciPy provides numerically stable evaluation of the CDF via scipy.stats.rice.cdf.
SciPy parameter mapping#
SciPy’s Rice distribution is scipy.stats.rice with shape parameter \(b=\nu/\sigma\):
def rice_logpdf(r, nu, sigma):
'''Log-PDF of Rice(nu, sigma) for r in [0, inf).
Uses exponentially-scaled Bessel i0e for numerical stability.
'''
r = np.asarray(r, dtype=float)
nu = float(nu)
sigma = float(sigma)
if nu < 0 or sigma <= 0:
raise ValueError("nu must be >= 0 and sigma must be > 0")
logpdf = np.full_like(r, -np.inf, dtype=float)
mask = r > 0
if not np.any(mask):
return logpdf
z = (r[mask] * nu) / (sigma**2)
# I0(z) grows ~ exp(z)/sqrt(2*pi*z), so work with i0e(z) = exp(-z) I0(z)
log_i0 = np.log(special.i0e(z)) + z
logpdf[mask] = (
np.log(r[mask])
- 2.0 * np.log(sigma)
- (r[mask] ** 2 + nu**2) / (2.0 * sigma**2)
+ log_i0
)
return logpdf
def rice_pdf(r, nu, sigma):
return np.exp(rice_logpdf(r, nu, sigma))
def rice_rvs_numpy(nu, sigma, size=1, rng=None):
'''Sample from Rice(nu, sigma) using NumPy only.
Construction: R = sqrt(X^2 + Y^2) where X~N(nu, sigma^2), Y~N(0, sigma^2).
'''
if rng is None:
rng = np.random.default_rng()
nu = float(nu)
sigma = float(sigma)
if nu < 0 or sigma <= 0:
raise ValueError("nu must be >= 0 and sigma must be > 0")
x = rng.normal(loc=nu, scale=sigma, size=size)
y = rng.normal(loc=0.0, scale=sigma, size=size)
return np.sqrt(x * x + y * y)
def ecdf(samples):
x = np.sort(np.asarray(samples, dtype=float))
y = np.arange(1, x.size + 1) / x.size
return x, y
4) Moments & Properties#
A convenient organizing principle is raw moments
Let \(a = \nu^2/(2\sigma^2)\).
Raw moments#
For \(n>-2\) (in particular for \(n=1,2,3,4\)):
where \({}_1F_1\) is the confluent hypergeometric function.
From these:
Mean: \(\mu = m_1\)
Variance: \(\mathrm{Var}(R)=m_2-m_1^2\)
A commonly used explicit expression for the mean (in Bessel functions) is
A very simple second raw moment comes from the 2D Gaussian construction:
Therefore:
Skewness and kurtosis#
From raw moments \(m_1,\dots,m_4\):
[ \mu_2 = m_2 - m_1^2,\quad \mu_3 = m_3 - 3 m_1 m_2 + 2 m_1^3,\quad \mu_4 = m_4 - 4 m_1 m_3 + 6 m_1^2 m_2 - 3 m_1^4. ]
[ \text{skewness} = \frac{\mu_3}{\mu_2^{3/2}},\qquad \text{kurtosis} = \frac{\mu_4}{\mu_2^{2}},\qquad \text{excess kurtosis} = \frac{\mu_4}{\mu_2^{2}} - 3. ]
SciPy can compute \((\mu,\mathrm{Var},\text{skew},\text{kurtosis})\) via stats.rice.stats(..., moments="mvsk").
MGF / characteristic function#
The MGF of \(R\) exists for all real \(t\) (the tail is essentially Gaussian), but it does not simplify to elementary functions:
A useful closed form exists for the squared amplitude \(R^2\):
\(R^2 = X^2+Y^2\) with \(X\sim\mathcal{N}(\nu,\sigma^2)\), \(Y\sim\mathcal{N}(0,\sigma^2)\).
Equivalently, \(R^2/\sigma^2 \sim \chi'^2(2,\lambda)\) with \(\lambda=(\nu/\sigma)^2\).
For \(t<1/(2\sigma^2)\):
The characteristic function follows by \(t\mapsto it\).
Entropy#
The differential entropy is
There is no simple closed form in general; it is commonly evaluated numerically (quadrature or Monte Carlo).
def rice_raw_moment(n, nu, sigma):
'''Raw moment m_n = E[R^n] for Rice(nu, sigma), using 1F1 representation.'''
nu = float(nu)
sigma = float(sigma)
if nu < 0 or sigma <= 0:
raise ValueError("nu must be >= 0 and sigma must be > 0")
a = nu**2 / (2.0 * sigma**2)
return (sigma**n) * (2.0 ** (n / 2.0)) * special.gamma(1.0 + n / 2.0) * special.hyp1f1(
-n / 2.0, 1.0, -a
)
def rice_mean(nu, sigma):
'''Mean of Rice(nu, sigma), computed in a numerically stable Bessel form.'''
nu = float(nu)
sigma = float(sigma)
if nu < 0 or sigma <= 0:
raise ValueError("nu must be >= 0 and sigma must be > 0")
a = nu**2 / (2.0 * sigma**2)
z = a / 2.0
# mean = sigma*sqrt(pi/2)*exp(-z)*[(1+a)I0(z) + a I1(z)]
# use i0e(z)=exp(-z)I0(z), i1e(z)=exp(-z)I1(z)
return sigma * np.sqrt(np.pi / 2.0) * ((1.0 + a) * special.i0e(z) + a * special.i1e(z))
def rice_central_moments(nu, sigma):
m1 = rice_raw_moment(1, nu, sigma)
m2 = rice_raw_moment(2, nu, sigma)
m3 = rice_raw_moment(3, nu, sigma)
m4 = rice_raw_moment(4, nu, sigma)
mu2 = m2 - m1**2
mu3 = m3 - 3 * m1 * m2 + 2 * m1**3
mu4 = m4 - 4 * m1 * m3 + 6 * m1**2 * m2 - 3 * m1**4
return m1, mu2, mu3, mu4
def rice_skew_kurtosis(nu, sigma):
mean, mu2, mu3, mu4 = rice_central_moments(nu, sigma)
skew = mu3 / (mu2 ** 1.5)
kurt = mu4 / (mu2**2)
return mean, mu2, skew, kurt
# Quick numeric check vs SciPy
nu0, sigma0 = 2.0, 1.3
b0 = nu0 / sigma0
mean_ours, var_ours, skew_ours, kurt_ours = rice_skew_kurtosis(nu0, sigma0)
mean_scipy, var_scipy, skew_scipy, kurt_scipy = stats.rice(b=b0, scale=sigma0).stats(
moments="mvsk"
)
(mean_ours, mean_scipy), (var_ours, var_scipy), (skew_ours, skew_scipy), (kurt_ours, kurt_scipy)
((2.474457072384077, 2.4744570723840766),
(1.2570621969284224, 1.2570621969284217),
(0.34234947509082275, 0.3423494750908338),
(2.835812529966139, -0.1641874700338617))
5) Parameter Interpretation#
Meaning of the parameters#
\(\nu\) (noncentrality / deterministic component)
Geometrically: distance of the 2D Gaussian mean from the origin.
In communications: proportional to the strength of the line-of-sight component.
\(\sigma\) (scale / scatter per quadrature)
Standard deviation of each Gaussian component \(X\) and \(Y\).
Controls how “spread out” the cloud is and therefore how variable the amplitude is.
A common dimensionless summary is the \(K\)-factor
interpretable as a (signal power)/(scatter power) ratio.
Shape changes#
Fix \(\sigma\) and increase \(\nu\):
density shifts to the right, mass near \(0\) decreases,
skewness decreases, and the distribution becomes more “Gaussian-like”.
Fix \(\nu\) and increase \(\sigma\):
density spreads out,
the mode moves left and the distribution becomes more right-skewed.
# How the PDF changes with nu (sigma fixed)
sigma_viz = 1.0
nus = [0.0, 0.5, 1.5, 3.0]
x = np.linspace(0, 8, 600)
fig = go.Figure()
for nu in nus:
fig.add_trace(
go.Scatter(x=x, y=rice_pdf(x, nu=nu, sigma=sigma_viz), mode="lines", name=f"ν={nu:g}")
)
fig.update_layout(
title=f"Rice PDF for varying ν (σ={sigma_viz:g} fixed)",
xaxis_title="r",
yaxis_title="density",
)
fig.show()
6) Derivations#
PDF (polar change of variables)#
Start from a shifted isotropic Gaussian:
The joint density is
Transform to polar coordinates \(x=r\cos\theta\), \(y=r\sin\theta\) with Jacobian \(|J|=r\):
Now integrate out the angle:
This yields the Rice PDF.
Expectation and variance#
The squared radius has a simple form:
Using \(X=\nu+\sigma Z_1\), \(Y=\sigma Z_2\) with \(Z_1,Z_2\sim\mathcal{N}(0,1)\) independent:
Then
Likelihood (i.i.d. sample)#
For data \(r_1,\dots,r_n\) (all \(\ge 0\)), the log-likelihood is
There is no closed-form MLE for \((\nu,\sigma)\) in general; numerical optimization is standard.
def rice_negloglik(params, data):
nu, sigma = float(params[0]), float(params[1])
if nu < 0 or sigma <= 0:
return np.inf
return -np.sum(rice_logpdf(data, nu, sigma))
# Synthetic data + MLE (numerical)
nu_true, sigma_true = 2.0, 1.0
n = 300
data = rice_rvs_numpy(nu_true, sigma_true, size=n, rng=rng)
# Start from SciPy fit as a reasonable initial guess (loc fixed at 0)
b_hat, loc_hat, scale_hat = stats.rice.fit(data, floc=0)
nu_init, sigma_init = b_hat * scale_hat, scale_hat
res = minimize(
rice_negloglik,
x0=np.array([nu_init, sigma_init]),
args=(data,),
method="L-BFGS-B",
bounds=[(0.0, None), (1e-12, None)],
)
nu_mle, sigma_mle = res.x
(nu_true, sigma_true), (nu_mle, sigma_mle), res.success
((2.0, 1.0), (1.9683706333582578, 0.9633278327527912), True)
7) Sampling & Simulation#
NumPy-only algorithm#
Use the defining construction:
Draw \(X\sim\mathcal{N}(\nu,\sigma^2)\).
Draw \(Y\sim\mathcal{N}(0,\sigma^2)\) independent.
Return \(R = \sqrt{X^2+Y^2}\).
Why this works: \(R\) is the radius of a 2D Gaussian centered at \((\nu,0)\).
This approach is fast, vectorized, and avoids special-function inversion.
# Smoke test: Monte Carlo mean/variance vs theory
nu_test, sigma_test = 1.5, 0.8
s = rice_rvs_numpy(nu_test, sigma_test, size=200_000, rng=rng)
mean_mc = s.mean()
var_mc = s.var()
mean_th = rice_mean(nu_test, sigma_test)
var_th = rice_raw_moment(2, nu_test, sigma_test) - mean_th**2
(mean_mc, mean_th), (var_mc, var_th)
((1.7337792125752307, 1.7345364036221205),
(0.5240102238815233, 0.5213834645096407))
8) Visualization#
We’ll compare:
PDF vs Monte Carlo histogram
CDF vs empirical CDF
We’ll compute the theoretical CDF using scipy.stats.rice (stable implementation).
nu_viz, sigma_viz = 2.0, 1.0
b_viz = nu_viz / sigma_viz
dist = stats.rice(b=b_viz, loc=0, scale=sigma_viz)
samples = rice_rvs_numpy(nu_viz, sigma_viz, size=80_000, rng=rng)
x_max = float(np.quantile(samples, 0.995))
x_grid = np.linspace(0, x_max, 600)
# PDF + histogram
fig = go.Figure()
fig.add_trace(
go.Histogram(
x=samples,
nbinsx=80,
histnorm="probability density",
name="Monte Carlo samples",
opacity=0.55,
)
)
fig.add_trace(
go.Scatter(
x=x_grid,
y=rice_pdf(x_grid, nu=nu_viz, sigma=sigma_viz),
mode="lines",
name="Theoretical PDF (NumPy+SciPy special)",
line=dict(width=3),
)
)
fig.update_layout(
title=f"Rice(ν={nu_viz:g}, σ={sigma_viz:g}): histogram vs PDF",
xaxis_title="r",
yaxis_title="density",
bargap=0.02,
)
fig.show()
# CDF + empirical CDF
xs, ys = ecdf(samples)
fig = go.Figure()
fig.add_trace(
go.Scatter(
x=x_grid,
y=dist.cdf(x_grid),
mode="lines",
name="Theoretical CDF (SciPy)",
line=dict(width=3),
)
)
fig.add_trace(
go.Scatter(
x=xs[::120],
y=ys[::120],
mode="markers",
name="Empirical CDF (subsampled)",
marker=dict(size=5),
)
)
fig.update_layout(
title=f"Rice(ν={nu_viz:g}, σ={sigma_viz:g}): empirical CDF vs CDF",
xaxis_title="r",
yaxis_title="CDF",
)
fig.show()
9) SciPy Integration (scipy.stats.rice)#
Key points:
Shape parameter:
b = nu / sigma(dimensionless).Scale parameter:
scale = sigma.Location shift:
loc(useloc=0for the standard Rice support \([0,\infty)\)).
Common methods:
pdf,logpdf,cdf,ppfrvs(random sampling)fit(MLE; consider fixingloc=0)
nu_true, sigma_true = 2.0, 1.0
b_true = nu_true / sigma_true
dist = stats.rice(b=b_true, loc=0, scale=sigma_true)
x = np.linspace(0, 8, 6)
pdf = dist.pdf(x)
cdf = dist.cdf(x)
samples_scipy = dist.rvs(size=5000, random_state=rng)
# Fit (fix loc=0 to match the usual Rice support)
b_hat, loc_hat, scale_hat = stats.rice.fit(samples_scipy, floc=0)
nu_hat = b_hat * scale_hat
(pdf, cdf), (nu_hat, scale_hat)
((array([0. , 0.346012, 0.250942, 0.012433, 0.000045, 0. ]),
array([0. , 0.242633, 0.841085, 0.995869, 0.99999 , 1. ])),
(1.9749664963783742, 0.9964787814644809))
10) Statistical Use Cases#
A) Hypothesis testing#
A common test in communications / detection:
\(H_0\): Rayleigh fading (\(\nu=0\))
\(H_1\): Rician fading (\(\nu>0\))
A likelihood ratio test compares the best fit under \(\nu=0\) (Rayleigh) versus the best fit under \(\nu\ge 0\) (Rice).
B) Bayesian modeling#
Rice likelihoods appear when you observe a magnitude but the latent signal is complex Gaussian. Typical Bayesian ingredients:
priors on \((\nu,\sigma)\) (or on \(K\) and \(\sigma\))
hierarchical models: different groups/locations share hyperparameters
inference via MCMC/VI since the posterior is non-conjugate
C) Generative modeling#
Rice is a natural emission distribution for nonnegative amplitudes generated from an underlying complex Gaussian latent:
sample latent complex value (mean + Gaussian noise)
emit its magnitude
This shows up in fading simulators and in synthetic MRI magnitude data generation.
# Simple likelihood-ratio statistic: Rayleigh (nu=0) vs Rice (nu>=0)
def rayleigh_sigma_mle(data):
data = np.asarray(data, dtype=float)
return np.sqrt(np.mean(data**2) / 2.0)
def loglik_rice(data, nu, sigma):
return float(np.sum(rice_logpdf(data, nu, sigma)))
def loglik_rayleigh(data, sigma):
return loglik_rice(data, nu=0.0, sigma=sigma)
# Simulate a dataset under H1
nu_alt, sigma_alt = 1.6, 1.0
x = rice_rvs_numpy(nu_alt, sigma_alt, size=250, rng=rng)
# Fit under H0 (Rayleigh) and H1 (Rice)
sigma0 = rayleigh_sigma_mle(x)
ll0 = loglik_rayleigh(x, sigma0)
b1, _, s1 = stats.rice.fit(x, floc=0)
nu1 = b1 * s1
ll1 = loglik_rice(x, nu1, s1)
llr = 2.0 * (ll1 - ll0)
(sigma0, (nu1, s1), llr)
(1.4926947250781553,
(1.5823631882104543, 0.9880517248571086),
12.58952187373336)
11) Pitfalls#
Invalid parameters:
Require \(\nu\ge 0\) and \(\sigma>0\) (SciPy:
b>=0,scale>0).
Parameterization mismatch:
Many references use \((\nu,\sigma)\); SciPy uses
b=nu/sigmaplusscale=sigma.Some domains use \(K=\nu^2/(2\sigma^2)\); make sure you track units.
locin SciPy:stats.rice.fitmay introduce a nonzerolocunless constrained.For the canonical Rice model, it’s common to fix
loc=0.
Numerical overflow in \(I_0\):
The Bessel function \(I_0(z)\) grows like \(\exp(z)\).
Prefer log-space computations and the scaled function
scipy.special.i0e.
Flat/ill-conditioned likelihoods:
When \(\nu/\sigma\) is very small, Rice is close to Rayleigh.
When \(\nu/\sigma\) is very large, \(R\) is close to normal around \(\nu\).
Both regimes can make numerical fitting sensitive to initialization.
# Numerical stability demo: naive PDF can overflow for large nu/sigma
nu_big, sigma_small = 30.0, 1.0
r_grid = np.linspace(0, 50, 6)
# Stable logpdf
lp = rice_logpdf(r_grid, nu_big, sigma_small)
# Naive pdf (may overflow internally via I0)
# Wrap in errstate: the point is to show unstable behavior without noisy warnings.
# NOTE: this is for demonstration; prefer rice_pdf / rice_logpdf above.
with np.errstate(over="ignore", under="ignore", invalid="ignore"):
naive = (
(r_grid / sigma_small**2)
* np.exp(-(r_grid**2 + nu_big**2) / (2 * sigma_small**2))
* special.i0((r_grid * nu_big) / (sigma_small**2))
)
lp, naive
(array([ -inf, -201.467827, -51.121463, -0.9188 , -50.774993,
-200.663442]),
array([ 0., 0., 0., nan, nan, nan]))
12) Summary#
Rice is a continuous distribution on \([0,\infty)\) modeling the magnitude of a shifted 2D Gaussian.
Parameters: \(\nu\) (deterministic/LOS component) and \(\sigma\) (scatter scale); \(K=\nu^2/(2\sigma^2)\) is a common dimensionless summary.
PDF involves a modified Bessel function \(I_0\); CDF can be written using the Marcum \(Q\)-function.
Raw moments have a compact form via \({}_1F_1\); \(\mathbb{E}[R^2]=2\sigma^2+\nu^2\) gives a simple variance identity.
Sampling is easy with NumPy: draw \((X,Y)\) as Gaussians and take \(\sqrt{X^2+Y^2}\).
References#
SciPy documentation:
scipy.stats.riceM. K. Simon and M.-S. Alouini, Digital Communication over Fading Channels (Rician fading)
Classic special-function relations for Rice / noncentral-\(\chi\) distributions